close all
clear all

yytot=0;
g=2.00E+16;
G=1.05E+08;
To=300;
Ceo=70*To;

%%%%% Read in data %%%%%

for nn=1:1
    fileloader1=sprintf('ratf%d=load(''sphere_nm2\\ratf%d.txt'');',nn,nn); % r a t f
    eval(fileloader1);
    fileloader2=sprintf('Tmax%d=load(''sphere_nm2\\ratf%d_max.txt'');',nn,nn); % Tmax
    eval(fileloader2);
    sizefind=sprintf('[yy%d,dummy]=size(ratf%d);',nn,nn);
    eval(sizefind);
    eval(sprintf('yytot=yytot+yy%d;',nn));
end

ratf=zeros(yytot,10);
yy0=0;
low=0;
high=0;

for nn=1:1
    low1=sprintf('low=1+high;');
    eval(low1);
    high1=sprintf('high=yy%d+low-1;',nn);
    eval(high1);
    mfiller1=sprintf('ratf(low:high,1:4)=ratf%d;',nn);
    eval(mfiller1);
    mfiller2=sprintf('ratf(low:high,8:10)=Tmax%d;',nn);
    eval(mfiller2);
end

%%%%% sort & arrange data %%%%%

epsilon=3*G./(g.*ratf(:,1));
beta=ratf(:,2).*ratf(:,4)./(Ceo*To*(4/3)*pi.*ratf(:,1).^3);
gamma=ratf(:,3)*g/Ceo;

ratf(:,5)=epsilon;
ratf(:,6)=beta;
ratf(:,7)=gamma;

temp1=sortrows(ratf,7);
temp2=sortrows(temp1,6);
ratfs=sortrows(temp2,5);

Tce=(ratfs(:,6)+1)*To;
Tcp=(ratfs(:,6)./ratfs(:,7)+1)*To;
Tcw=(ratfs(:,6)./(ratfs(:,7).*ratfs(:,5))+1)*To;

thetae=(ratfs(:,8)-To)./(Tce-To);
thetap=(ratfs(:,9)-To)./(Tcp-To);
thetaw=(ratfs(:,10)-To)./(Tcw-To);

    count=1;
    deps=diff(ratfs(1:yytot,5));
    for tt=1:yytot-1
        if deps(tt)>0
            count=count+1;
        end
    end
    neps=count;
    leps=yytot/neps;
    
    count=1;
    dbeta=diff(ratfs(2*leps+1:3*leps,6));
    for tt=1:leps-1
        if dbeta(tt)>0
            count=count+1;
        end
    end
    nbeta=count;
    lbeta=leps/nbeta;
    ngamma=lbeta;
    
    vgamma=gamma(1:ngamma);
    vbeta=zeros(nbeta,1);
    for kk=1:nbeta
        vbeta(kk,1)=beta(kk*lbeta-1,1);
    end
    
        veps=zeros(neps,1);
    for kk=1:neps
        veps(kk,1)=epsilon(kk*yytot/neps-1,1);
    end
    count=0;
    for kk=1:neps
        eval(sprintf('T%d=zeros(ngamma,nbeta,3);',kk));
        for oo=1:nbeta
            count=count+1;
            olow=(count-1)*ngamma+1;
            ohigh=(count)*ngamma;
            eval(sprintf('T%d(1:ngamma,oo,1)=thetae(olow:ohigh,1);',kk));
            eval(sprintf('T%d(1:ngamma,oo,2)=thetap(olow:ohigh,1);',kk));
            eval(sprintf('T%d(1:ngamma,oo,3)=thetaw(olow:ohigh,1);',kk));
        end
    end
    
%%%%% Creating cut-off vectors for fitting %%%%%

    % observed fit cutoffs
    %Te: gamma < 1000
    %Tp: gamma < 300
    %Tw: gamma < 1000
    egc=1000;
    pgc=300;
    wgc=1000;
    
    [yegct,dummy]=find(vgamma<=egc);yegc=length(yegct);
    [ypgct,dummy]=find(vgamma<=pgc);ypgc=length(ypgct);
    [ywgct,dummy]=find(vgamma<=wgc);ywgc=length(ywgct);
    segamma=vgamma(1:yegc);
    spgamma=vgamma(1:ypgc);
    swgamma=vgamma(1:ywgc);
    
    for kk=1:neps
        eval(sprintf('Tegc%d=T%d(1:yegc,:,1);',kk,kk));
        eval(sprintf('Tpgc%d=T%d(1:ypgc,:,2);',kk,kk));
        eval(sprintf('Twgc%d=T%d(1:ywgc,:,3);',kk,kk));
    end

    sle=neps*nbeta*yegc;
    slp=neps*nbeta*ypgc;
    slw=neps*nbeta*ywgc;
    
    gse=zeros(sle,1);bse=zeros(sle,1);ese=zeros(sle,1);
    gsp=zeros(slp,1);bsp=zeros(slp,1);esp=zeros(slp,1);
    gsw=zeros(slw,1);bsw=zeros(slw,1);esw=zeros(slw,1);
    
    thetae_s=gse;thetap_s=gsp;thetaw_s=gsw;
    count=0;

    for kk=1:neps
        for oo=1:nbeta
        %oo=1;
            count=count+1;
            oloe=(count-1)*yegc+1;ohie=(count)*yegc;
            olop=(count-1)*ypgc+1;ohip=(count)*ypgc;
            olow=(count-1)*ywgc+1;ohiw=(count)*ywgc;
            
            btemp=vbeta(oo,1);
            etemp=veps(kk,1);
            gse(oloe:ohie,1)=vgamma(1:yegc,1);bse(oloe:ohie,1)=btemp;ese(oloe:ohie,1)=etemp;
            gsp(olop:ohip,1)=vgamma(1:ypgc,1);bsp(olop:ohip,1)=btemp;esp(olop:ohip,1)=etemp;
            gsw(olow:ohiw,1)=vgamma(1:ywgc,1);bsw(olow:ohiw,1)=btemp;esw(olow:ohiw,1)=etemp;
            
            eval(sprintf('thetae_s(oloe:ohie,1)=Tegc%d(1:yegc,oo);',kk));
            eval(sprintf('thetap_s(olop:ohip,1)=Tpgc%d(1:ypgc,oo);',kk));
            eval(sprintf('thetaw_s(olow:ohiw,1)=Twgc%d(1:ywgc,oo);',kk));
            
        end

%         scatter3(gse(1:ohie),bse(1:ohie),thetae_s(1:ohie));
%         set(gca, 'xscale', 'log','yscale', 'log','zscale', 'log')
%         xlabel('\gamma');
%         ylabel('\beta');
%         hold on
    end
%    scatter3(gse(1:yegc),bse(1:yegc),thetae_s(1:yegc));
%    set(gca, 'xscale', 'log','yscale', 'log','zscale', 'log')
%    xlabel('\gamma');
%    ylabel('\beta');
%    hold on
   %scatter3(gse(yegc+1:2*yegc),bse(yegc+1:2*yegc),thetae_s(yegc+1:2*yegc),'r*');
   %scatter3(gse,bse,thetae_s);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%% Curve fitting %%%%

logte=log(thetae_s);
epara(:,3)=log(ese);
epara(:,2)=log(bse);
epara(:,1)=log(gse);

logtp=log(thetap_s);
ppara(:,3)=log(esp);
ppara(:,2)=log(bsp);
ppara(:,1)=log(gsp);

logtw=log(thetaw_s);
wpara(:,3)=log(esw);
wpara(:,2)=log(bsw);
wpara(:,1)=log(gsw);

options.MaxFunEvals=1e8;
options.TolX = 1e-7;
options.TolFun = 1e-15;

%%%% lattice C0 %%%%
CP=[-7.5;1;10000;9;1;-0.2;-0.15;1.6];
%%%% water C0 %%%%
CW=[-8.4;0.6;1;10000;9;1;-0.2;-0.12;1.6;0.001];
%%%% electron C0 %%%%
CE=[-2.1;0.7;0.7;0;-0.25;-0.35;0.7;1.4;0;1;-0.13;-0.07;-0.11;1;-3;pi/2;1;-3;pi/2];

[c,res]=lsqcurvefit(@Tfunce,CE,epara,logte);
cfite=c;
[c,res]=lsqcurvefit(@Tfuncp,CP,ppara,logtp);
cfitp=c;
[c,res]=lsqcurvefit(@Tfuncw,CW,wpara,logtw);
cfitw=c;

logte_out=Tfunce(cfite,epara);
logtp_out=Tfuncp(cfitp,ppara);
logtw_out=Tfuncw(cfitw,wpara);

res_ec=0;
res_pc=0;
res_wc=0;

for ee=1:sle
    res_ec=((logte(ee)-logte_out(ee))/logte(ee))^2+res_ec;
end
for ee=1:slp
    res_pc=((logtp(ee)-logtp_out(ee))/logtp(ee))^2+res_pc;
end
for ee=1:slw
    res_wc=((logtw(ee)-logtw_out(ee))/logtw(ee))^2+res_wc;
end

res_e=(res_ec/sle)^(1/2);
res_p=(res_pc/slp)^(1/2);
res_w=(res_wc/slw)^(1/2);

figure(1)
scatter3(gse,bse,exp(logte_out));
set(gca, 'xscale', 'log','yscale', 'log','zscale', 'log')
xlabel('\gamma');
ylabel('\beta');
hold on
scatter3(gse,bse,thetae_s);

figure(2)
scatter3(gsp,bsp,exp(logtp_out));
set(gca, 'xscale', 'log','yscale', 'log','zscale', 'log')
xlabel('\gamma');
ylabel('\beta');
hold on
scatter3(gsp,bsp,thetap_s);

figure(3)
scatter3(gsw,bsw,exp(logtw_out));
set(gca, 'xscale', 'log','yscale', 'log','zscale', 'log')
xlabel('\gamma');
ylabel('\beta');
hold on
scatter3(gsw,bsw,thetaw_s);

thetaw_max=max(thetaw);
thetap_max=max(thetap);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
pmode=7;


if pmode==4;
    pl=lbeta;
    close all
    %%% SCATTER PLOTS %%%
for ss=1:yytot/pl-1
    plow=pl*(ss-1)+1;
    phigh=pl*ss;
    figure(3);
    xlabel('\Gamma','fontsize',14);
    ylabel('\beta','fontsize',14);
    zlabel('\theta_w','fontsize',14);
    if ratfs(plow,5)==ratfs(leps*0+1,5);
        scatter3(ratfs(plow:phigh,7),ratfs(plow:phigh,6),thetaw(plow:phigh,1),'r.');
    elseif ratfs(plow,5)==ratfs(leps*1+1,5);
        scatter3(ratfs(plow:phigh,7),ratfs(plow:phigh,6),thetaw(plow:phigh,1),'b.');
    elseif ratfs(plow,5)==ratfs(leps*2+1,5);
        scatter3(ratfs(plow:phigh,7),ratfs(plow:phigh,6),thetaw(plow:phigh,1),'g.');
    elseif ratfs(plow,5)==ratfs(leps*3+1,5);
        scatter3(ratfs(plow:phigh,7),ratfs(plow:phigh,6),thetaw(plow:phigh,1),'c.');
    elseif ratfs(plow,5)==ratfs(leps*4+1,5);
        scatter3(ratfs(plow:phigh,7),ratfs(plow:phigh,6),thetaw(plow:phigh,1),'m.');
    end
    set(gca, 'xscale', 'log','yscale', 'log','zscale', 'log')
    axis tight
    %set(gca, 'zscale', 'log')
    hold on
end

for ss=1:yytot/pl-1
    plow=pl*(ss-1)+1;
    phigh=pl*ss;
    figure(4);
    xlabel('\Gamma','fontsize',14);
    ylabel('\beta','fontsize',14);
    zlabel('\theta_p','fontsize',14);
    if ratfs(plow,5)==ratfs(leps*0+1,5);
        scatter3(ratfs(plow:phigh,7),ratfs(plow:phigh,6),thetap(plow:phigh,1),'r.');
    elseif ratfs(plow,5)==ratfs(leps*1+1,5);
        scatter3(ratfs(plow:phigh,7),ratfs(plow:phigh,6),thetap(plow:phigh,1),'b.');
    elseif ratfs(plow,5)==ratfs(leps*2+1,5);
        scatter3(ratfs(plow:phigh,7),ratfs(plow:phigh,6),thetap(plow:phigh,1),'g.');
    elseif ratfs(plow,5)==ratfs(leps*3+1,5);
        scatter3(ratfs(plow:phigh,7),ratfs(plow:phigh,6),thetap(plow:phigh,1),'c.');
    elseif ratfs(plow,5)==ratfs(leps*4+1,5);
        scatter3(ratfs(plow:phigh,7),ratfs(plow:phigh,6),thetap(plow:phigh,1),'m.');
    end
    set(gca, 'xscale', 'log','yscale', 'log','zscale', 'log')
    axis tight
    %set(gca, 'zscale', 'log')
    hold on
end

for ss=1:yytot/pl-1
    plow=pl*(ss-1)+1;
    phigh=pl*ss;
    figure(5);
    xlabel('\Gamma','fontsize',14);
    ylabel('\beta','fontsize',14);
    zlabel('\theta_e','fontsize',14);
    if ratfs(plow,5)==ratfs(leps*0+1,5);
        scatter3(ratfs(plow:phigh,7),ratfs(plow:phigh,6),thetae(plow:phigh,1),'r.');
    elseif ratfs(plow,5)==ratfs(leps*1+1,5);
        scatter3(ratfs(plow:phigh,7),ratfs(plow:phigh,6),thetae(plow:phigh,1),'b.');
    elseif ratfs(plow,5)==ratfs(leps*2+1,5);
        scatter3(ratfs(plow:phigh,7),ratfs(plow:phigh,6),thetae(plow:phigh,1),'g.');
    elseif ratfs(plow,5)==ratfs(leps*3+1,5);
        scatter3(ratfs(plow:phigh,7),ratfs(plow:phigh,6),thetae(plow:phigh,1),'c.');
    elseif ratfs(plow,5)==ratfs(leps*4+1,5);
        scatter3(ratfs(plow:phigh,7),ratfs(plow:phigh,6),thetae(plow:phigh,1),'m.');
    end
    set(gca, 'xscale', 'log','yscale', 'log','zscale', 'log')
    axis tight
    %set(gca, 'zscale', 'log')
    hold on
end
% water
% scatter3(1300,1,0.6,'mo');
% scatter3(6000,1,1.5,'co');
% scatter3(30000,1,3,'go');
% phonon
%  scatter3(1788,1,0.5,'mo');
%  scatter3(5059,1,1.5,'co');
%  scatter3(20000,1,6,'go');
%  scatter3(56563,1,20,'bo');

%2*gamma*epsilon pplateu bounds

%%% CONTOUR PLOTS %%%
elseif pmode==5;
    Tmatrix=zeros(lbeta,nbeta);
    betavector=zeros(nbeta,1);
    gammavector=ratfs(1:lbeta,7);
    for tt=1:nbeta
        low=(tt-1)*lbeta+1;
        high=tt*lbeta;
        betavector(tt,1)=ratfs(low,6);
        Tmatrix(:,tt)=thetae(low:high,1);
    end
    
    figure(3);
    xlabel('\Gamma');
    ylabel('\beta');
    zlabel('\theta_e');
    %surf(Tmatrix);
    surf(betavector,gammavector,Tmatrix);
    set(gca, 'xscale', 'log','yscale', 'log','zscale', 'log')
%     elseif ratfs(plow,5)==ratfs(floor(yytot/2),5);
%         loglog(ratfs(plow:phigh,7),thetaw(plow:phigh,1),'bo');
%     elseif ratfs(plow,5)==ratfs(yytot,5);
%         loglog(ratfs(plow:phigh,7),thetaw(plow:phigh,1),'go');
    hold on
    %surfc(ratfs(eta_lim+1:2*eta_lim,7),ratfs(eta_lim+1:2*eta_lim,6),thetaw(eta_lim+1:2*eta_lim,1));
    %surfc(ratfs(eta_lim+1:2*eta_lim,7),ratfs(eta_lim+1:2*eta_lim,6),thetaw(eta_lim+1:2*eta_lim,1));
end
%save('sphere\\ratf_tot.txt','ratfs','-ASCII','-TABS');
mode=1;
if mode==2
    para(:,3)=log(epsilon);
    para(:,2)=log(beta);
    para(:,1)=log(gamma);
    % para(:,3)=epsilon;
    % para(:,2)=beta;
    % para(:,1)=gamma;
    log_tp=log(thetap);
    log_tw=log(thetaw);
    log_te=log(thetae);
    options.MaxFunEvals=1e8;
    options.TolX = 1e-7;
    options.TolFun = 1e-15;

    % lattice C0 %
    C0=[-7.5;1;10000;9;1;-0.2;-0.15;1.6];
    % water C0 %
    %C0=[-8.4;0.6;1;10000;9;1;-0.2;-0.12;1.6];  %lsq able
    %C0=[-7;0.6;0.75;10000;9;1;-0.2;-0.12;1.6;1;100;0.5;-2;7;0.2];
    % electron C0 %
    %C0=[-2.1;0.7;0.7;0;-0.25;-0.35;0.7;1.4;0;1;-0.13;-0.07;-0.11;1;-3;pi/2;1;-3;pi/2];
    %[c,res]=lsqcurvefit(@Tfunce,C0,para,log_te);
    log_t_out=Tfuncp(C0,para);

    res0t=0;
    for ee=1:yytot
        res0t=(log_tw(ee)-log_t_out(ee))^2+res0t;
    end
    res0=(res0t/yytot)^(1/2);

    if pmode==4
        %scatter3(gamma,beta,exp(log_t_out),'k.');
    end
end
% delta=abs(log_te-log_t_out);
% figure(2)
% scatter3(gamma,beta,exp(delta),'b.');
% set(gca, 'xscale', 'log','yscale', 'log');%,'zscale', 'log')

iflag=0;
if iflag==1
%%% ITERATION FOR BOIL/MELT/EXPLODE %%%
Tb=647;
Tm=1337;

betam=zeros(lbeta,1);
betab=zeros(lbeta,1);
betae=zeros(lbeta,1);
test_gamma=gamma(1:lbeta,1);
test_eps=3*G/(25e-9*g);

delta_tol1=0.2;
delta_tol2=0.005;
step=(beta(lbeta+1,1)-beta(1,1))*100;

% water C0 %
C0=[-8.3;0.75;1;10000;9;1;-0.2;-0.15;1.6];

for tt=1:lbeta
    delta=1;
    itno=0;
    guess_beta=min(beta);
    idsp=sprintf('disp(''Now iterating on gamma %d'');',test_gamma(tt,1));
    eval(idsp);
    while delta>delta_tol1
        itno=itno+1;
        %figure(1)
        %loglog(itno,delta)
        %hold on
        
        guess_beta=guess_beta+step;
        
        Tcee=(guess_beta+1)*To;
        Tcpm=(guess_beta./test_gamma(tt)+1)*To;
        Tcwb=(guess_beta./(test_gamma(tt).*test_eps)+1)*To;

        Tee=(10000-To)/(Tcee-To);
        Tpm=(Tm-To)/(Tcpm-To);
        Twb=(Tb-To)/(Tcwb-To);

        ip(1,3)=log(test_eps);
        ip(1,2)=log(guess_beta);
        ip(1,1)=log(test_gamma(tt));

        log_t_out=Tfuncw(C0,ip);
        delta=abs((log(Twb)-log_t_out)/log(Twb));
    end
    while delta>delta_tol2
        itno=itno+1;
        %figure(1)
        %loglog(itno,delta)
        %hold on
        guess_beta=guess_beta+step/5000;
        
        Tcee=(guess_beta+1)*To;
        Tcpm=(guess_beta./test_gamma(tt)+1)*To;
        Tcwb=(guess_beta./(test_gamma(tt).*test_eps)+1)*To;

        Tee=(10000-To)/(Tcee-To);
        Tpm=(Tm-To)/(Tcpm-To);
        Twb=(Tb-To)/(Tcwb-To);

        ip(1,3)=log(test_eps);
        ip(1,2)=log(guess_beta);
        ip(1,1)=log(test_gamma(tt));

        log_t_out=Tfuncw(C0,ip);
        delta=abs((log(Twb)-log_t_out)/log(Twb));
    end
    betab(tt,1)=guess_beta;
end

delta_tol1=0.2;
delta_tol2=0.005;
step=(beta(lbeta+1,1)-beta(1,1))*100;

% lattice C0 %
C0=[-7.5;1;10000;9;1;-0.2;-0.15;1.6];

for tt=1:lbeta
    delta=1;
    itno=0;
    guess_beta=min(beta);
    idsp=sprintf('disp(''Now iterating on gamma %d'');',test_gamma(tt,1));
    eval(idsp);
    while delta>delta_tol1
        itno=itno+1;
        %figure(1)
        %loglog(itno,delta)
        %hold on
        
        guess_beta=guess_beta+step;
        
        Tcee=(guess_beta+1)*To;
        Tcpm=(guess_beta./test_gamma(tt)+1)*To;
        Tcwb=(guess_beta./(test_gamma(tt).*test_eps)+1)*To;

        Tee=(10000-To)/(Tcee-To);
        Tpm=(Tm-To)/(Tcpm-To);
        Twb=(Tb-To)/(Tcwb-To);

        ip(1,3)=log(test_eps);
        ip(1,2)=log(guess_beta);
        ip(1,1)=log(test_gamma(tt));

        log_t_out=Tfuncp(C0,ip);
        delta=abs((log(Tpm)-log_t_out)/log(Tpm));
    end
    while delta>delta_tol2
        itno=itno+1;
        %figure(1)
        %loglog(itno,delta)
        %hold on
        guess_beta=guess_beta+step/5000;
        
        Tcee=(guess_beta+1)*To;
        Tcpm=(guess_beta./test_gamma(tt)+1)*To;
        Tcwb=(guess_beta./(test_gamma(tt).*test_eps)+1)*To;

        Tee=(10000-To)/(Tcee-To);
        Tpm=(Tm-To)/(Tcpm-To);
        Twb=(Tb-To)/(Tcwb-To);

        ip(1,3)=log(test_eps);
        ip(1,2)=log(guess_beta);
        ip(1,1)=log(test_gamma(tt));

        log_t_out=Tfuncp(C0,ip);
        delta=abs((log(Tpm)-log_t_out)/log(Tpm));
    end
    betam(tt,1)=guess_beta;
end

delta_tol1=0.2;
delta_tol2=0.01;
step=(beta(lbeta+1,1)-beta(1,1))*100;

% electron C0 %
C0=[-2.1;0.7;0.7;0;-0.25;-0.35;0.7;1.4;0;1;-0.13;-0.07;-0.11;1;-3;pi/2;1;-3;pi/2];

for tt=1:lbeta
    delta=1;
    itno=0;
    guess_beta=min(beta);
    idsp=sprintf('disp(''Now iterating on gamma %d'');',test_gamma(tt,1));
    eval(idsp);
    while delta>delta_tol1
        itno=itno+1;
        %figure(1)
        %loglog(itno,delta)
        %hold on
        
        guess_beta=guess_beta+step;
        
        Tcee=(guess_beta+1)*To;
        Tcpm=(guess_beta./test_gamma(tt)+1)*To;
        Tcwb=(guess_beta./(test_gamma(tt).*test_eps)+1)*To;

        Tee=(10000-To)/(Tcee-To);
        Tpm=(Tm-To)/(Tcpm-To);
        Twb=(Tb-To)/(Tcwb-To);

        ip(1,3)=log(test_eps);
        ip(1,2)=log(guess_beta);
        ip(1,1)=log(test_gamma(tt));

        log_t_out=Tfunce(C0,ip);
        delta=abs((log(Tee)-log_t_out)/log(Tee));
    end
    while delta>delta_tol2
        itno=itno+1;
        %figure(1)
        %loglog(itno,delta)
        %hold on
        guess_beta=guess_beta+step/1000;
        
        Tcee=(guess_beta+1)*To;
        Tcpm=(guess_beta./test_gamma(tt)+1)*To;
        Tcwb=(guess_beta./(test_gamma(tt).*test_eps)+1)*To;

        Tee=(10000-To)/(Tcee-To);
        Tpm=(Tm-To)/(Tcpm-To);
        Twb=(Tb-To)/(Tcwb-To);

        ip(1,3)=log(test_eps);
        ip(1,2)=log(guess_beta);
        ip(1,1)=log(test_gamma(tt));

        log_t_out=Tfunce(C0,ip);
        delta=abs((log(Tee)-log_t_out)/log(Tee));
    end
    betae(tt,1)=guess_beta;
end

figure(11)
semilogx(test_gamma,betam,'r-');
hold on
semilogx(test_gamma,betab,'b-');
semilogx(test_gamma(1:17,1),betae(1:17,1),'g-');
xlabel('\Gamma');
ylabel('\beta');
legend('Particle Melt Initiation','Water Super Critical','Coulombic Explosion');
axis([10e-3 10e5 2500 10000]);
end